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In genetic circuits, when the mRNA lifetime is short compared to the cell cycle, proteins are 
produced in geometrically-distributed bursts, which greatly affects the cellular switching dynamics 
between different metastable phenotypic states. Motivated by this scenario, we study a general 
problem of switching or escape in stochastic populations, where influx of particles occurs in groups 
or bursts, sampled from an arbitrary distribution. The fact that the step size of the influx reac¬ 
tion is a-priori unknown, and in general, may fluctuate in time with a given correlation time and 
statistics, introduces an additional non-demographic step-size noise into the system. Employing the 
probability generating function technique in conjunction with Hamiltonian formulation, we are able 
to map the problem in the leading order onto solving a stationary Hamilton-Jacobi equation. We 
show that bursty influx exponentially decreases the mean escape time compared to the “usual case” 
of single-step influx. In particular, close to bifurcation we find a simple analytical expression for the 
mean escape time, which solely depends on the mean and variance of the burst-size distribution. 
Our results are demonstrated on several realistic distributions and compare well with numerical 
Monte-Carlo simulations. 
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Stochastic populations possessing multiple long-lived 
metastable states can undergo noise-induced switching, 
which plays a key role in many systems in physics, chem¬ 
istry, biology, and ecology [1-4]. Many previous studies 
of such noise-induced switching or escape have focused on 
the role of demographic noise, see e.g., Refs. [5-14]. Here 
a stochastic population, which regulates itself through 
random births and deaths, eventually escapes its initial 
metastable state via a rare sequence of events caused by 
the demographic noise. 

In models of stochastic population dynamics, an im¬ 
portant process that prevents the population from going 
extinct and is the main driver towards establishment of 
the population is the process of immigration or influx of 
particles [14-16]. In most cases, immigration is modeled 
using single-step influx (SSI). Here, influx occurs via a 
single-step process 0 —A, which corresponds to a sin¬ 
gle individual arriving at a certain rate, increasing the 
population size from n to n -I- 1 individuals. 

Yet, a more general and realistic way to model immi¬ 
gration is to consider bursty influx (HI) of particles or 
individuals. Here, influx occurs in bursts of individuals, 
with the burst size sampled from an arbitrary burst-size 
distribution (BSD). This generalized BI process appears 
in various scientific areas. In ecology, BI can account 
for more complex migration patterns such as arrival in 
groups [13, 17]. In gene expression, when the mRNA 
lifetime is short compared to the cell cycle, protein 
synthesis effectively occurs in geometrically-distributed 
bursts [18-27]. This gives rise to non-Poissonian pro¬ 
tein statistics [20, 22, 26, 27] and may exponentially 
decrease switching times between different phenotypic 
states [21, 23, 24, 27, 28]. Moreover, in evolutionary biol¬ 
ogy it has been shown that rapid punctuational bursts of 
speciation is a significant mechanism responsible for the 
evolution of species [29]. A similar mechanism of BI of 


lexical changes is also key in evolutionary linguistics [30] . 

Importantly, accounting for uncertainty in the reaction 
step size introduces additional non-demographic noise 
into the system, which can be called step-size noise, and 
has not been systematically investigated before. This 
noise, in general, includes a step-size parameter, k{t), 
fluctuating in time, with a given statistics and correla¬ 
tion time. Incorporating reaction step-size noise is ex¬ 
pected to dramatically change the statistics of interest. 
This is since, e.g., rare events such as switching or escape, 
are expected to be determined by a competition between 
the tails of the probability distribution of the population 
sizes given a certain burst size, and the tails of the BSD. 

In this paper, we make a first step towards a systematic 
investigation of generic reaction step-size noise. We do 
so by considering static step-size noise in the form of BI, 
0 —>■ kA, where k can take random integer values drawn 
from an arbitrary BSD. We demonstrate our formalism 
on a specific Verhulst-like model which includes, apart 
from BI, binary reproduction and death. This model 
gives rise to two fixed points at the deterministic level: a 
lower stable point ni, and a higher unstable point n 2 , see 
below. As a result, the stochastic population which ini¬ 
tially resides in the vicinity of the long-lived metastable 
state Til, ultimately crosses the unstable point n 2 and es¬ 
capes to infinity, leading to a Malthusian catastrophe [5] . 
In order to study how BI affects the mean escape time 
(MET) as function of the BSD, we employ the probabil¬ 
ity generating function formalism [3]. By doing so, we 
transform the master equation into a single partial dif¬ 
ferential equation. To this end, we use the WKB theory 
to treat this equation and find the MET in the leading 
exponential order as function of the BSD. Our analy¬ 
sis shows that BI can exponentially decrease the MET 
compared to the usual SSI case. In the bifurcation limit 
where the stable and unstable fixed points are close, we 
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derive a simple expression for the reduction of the MET, 
which depends solely on the first and second moments of 
the BSD. Our results are demonstrated on several BSD 
examples such as the geometric, negative-binomial, and 
scale-free distributions, all of which compare well with 
numerical Monte-Carlo simulations. 

Let n(t) be the population size. Our model includes bi¬ 
nary reproduction 2A 3 ^ death A ^ 0, where 
a > 1 is the reproduction rate, N :$> 1 is related to the 
critical population size, see below, and time is rescaled 
by the death rate. Accounting for immigration via a 

single-step influx (SSI) process 0 A, where ko is the 
immigration rate, we obtain the following deterministic 
rate equation for the mean population size n{t) 

- 7 

n= ko + - n. (I) 

In order to study bursty influx (BI) we replace the 
SSI process by an infinite set of influx processes 0 —>■ kA 
(k = 0, 1 , 2 ,....), where the burst size, k, is drawn from 
an arbitrary BSD, Dlk), which satisfies ~ 

and has a finite mean (A:) and variance varik). The rate 
of the BI processes 0 —>■ kA is taken to be kQD{k)/{k). 
This insures, see below, that the rate equation ( 1 ) and 
the corresponding fixed points are unaffected by the BI. 

Eq. (1) has two fixed points. The lower, ni/N =(1 — 
— 4:ako/N)/{2a), is an attracting fixed point, and the 
higher, n 2 /N = (1 -|- — 4:ako/N)/{2a), is a repelling 

fixed point, where 712 > ni, and we assume that ni 
1. This system describes a runaway model. Here an 
initial population, no < 712 , first relaxes into the long- 
lived metastable state, tii. Yet, after an exponentially 
long waiting time, the population crosses the barrier at 
the critical population size n 2 , and reaches a state of 
unlimited growth [31]. Our aim here is to investigate 
how the MET depends on the functional form of D{k). 

To do so, we account for demographic stochasticity and 
BI, and consider the corresponding master equation for 
Pn{t) - the probability to find n individuals at time t 


Multiplying Eq. (2) by p", and summing over all values of 
n we arrive at a single evolution equation for G(p, t) [35] 

^ - 1) b - 

+ ^p^[p-l)dlG. (3) 

Here we have introduced g{p) = This 

function serves as the burst-size probability generating 
function and encodes the arbitrary (known) BSD. 

As the MET is exponentially large, see below, once 
the system has settled into the vicinity of the long-lived 
metastable state tii we can replace dtG by zero in Eq. (3). 
The quasi-stationary solution of Eq. (3), Gqsip), with 
proper boundary conditions encodes the quasi-stationary 
distribution of the population sizes around ni, which 
yields the MET in the leading order [32-34], whereas 
subleading order calculations require spectral analysis of 
Eq. (3) [33] . The ordinary differential equation for Gqs (p) 
is, in general, unsolvable analytically, and we thus use the 
WKB theory in the spirit of [34]. Employing the WKB 
ansatz Gqs{p) = in Eq. (3), where S{p) is the ac¬ 

tion, we obtain in the leading 0{N) order, a stationary 
Hamilton-Jacobi equation, 'H[p, —S'{p)] = 0, with 

Pip. q) = {p- 1) [koUp) -q + ap^g^] . (4) 

Here we have defined kq = k^/N and /(p) = [p(p) — 
!)]• Note, that /(p) does not diverge at p = 1 
and satisfies /(I) = 1. In (4) we have further introduced 
g(p) = —S'{p), which plays the role of the reaction coor¬ 
dinate conjugate to the momentum p [32, 34]. 

The equation 'H = 0 has three solutions; one of them, 
p = 1, is the trivial solution called the relaxation trajec¬ 
tory, and corresponds to the deterministic dynamics [32] . 
Indeed, at p = I, the Hamilton equation for q satisfies 
q = Ko + aq^ — q, which coincides with rate equation (1), 
upon defining q = n/N as the population concentration. 
The Hamiltonian, Eq. (4), has two more, non-trivial. 


= 


(fc) 

-k 


.k=0 

a{n — 1)^ 


fe =0 

2 


( 2 ) 


an 


Pn-l -+ (77 -I- l)P„+i — nPm 


where the terms in square brackets correspond to the BI. 
Multiplying Eq. (2) by n and summing over all values 
of 71, we recover rate equation (1), justifying a-posteriori 
the rate chosen for the BI process. That Eq. (2) includes 
reactions with all possible step sizes renders this problem 
difficult to treat via the real-space approach [7-12]. We 
thus use the momentum-space approach [32-34] instead. 

We now define the probability generating function of 
the population, G = J2^=oP'^^ri [3], where p is an auxil¬ 
iary variable which plays the role of the momentum [32] . 



FIG. 1. (Color online) Shown are characteristic activation 
trajectories q+ (p) (solid line) and q- (p) (dashed line) given 
by Eq. (5), in the case of geometric BSD. The shaded area 
represents the accumulated action, see Eq. (6). Here param¬ 
eters are kq = 0.2, a = 0.3, and b — 3, see Table I. 
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zero-energy trajectories, called the activation trajectories 


<l±{p) 


l±^yl- 4aKop'^f{p) 
2ap^ 


(5) 


Figure 1 shows characteristic activation trajectories 
between the two fixed points, qi ^2 = ''^ 1 , 2 /-^, hr the case 
of geometric BSD. Having found the activation trajecto¬ 
ries, the MET r is given in the leading order by [32, 37] 


NAS 


AS = 


[<}+{p) - <l-ip)]dp- ( 6 ) 



b 1 


The accumulated action, AS, corresponds to the area 
between the activation and relaxation trajectories, see 
Fig. 1, and q+{p) and q-{p) are determined by Eq. (5). 
In addition, Pmax is determined by the condition 
g+(pmax) = 9 -(Pmax) which yields a transcendental eqna- 
tion Pmax/(Pmax) = 1/ (dcrfvQj. The condition for Pniax 
stems from the monotonicity in p of the trajectories. 
Eq. ( 6 ) for the MET together with Eq. (5) in the case 
of an arbitrary BSD is one of our main results. 

To evaluate the effect BI has on the MET we compare 
Eq. ( 6 ) to the simple case of SSI, see Table I. Here, AS 
from Eq. ( 6 ) can be explicitly computed and we find 


ASssi 



4-1-2 arcsin(2Y'aKo) — tJ" 


(7) 


Therefore, the general result for r for an arbitrary BSD 
can be written as r ~ = (tssi)^: where we define 

(j) = AS/ ASssi as the ratio between the accumulated 
actions in the BI and SSI cases. In the following we 
calculate the exponent (p for various BSDs. 

Figure 2 compares between the theoretical and numer¬ 
ical results for the MET, for different BSDs defined in 
Table I. Here, our numerical simulations were performed 
using an extended version of the Gillespie algorithm [38] 
that accounts for BI. For all BSDs considered, an expo¬ 
nential reduction of the MET is observed. To remind 



D(k) 

(k) 

var(fc) 

f(p) 

SSI 

^k,l 

1 

0 

1 

KSI 

Sk,K 

K 

0 

pK-i 

K(p-l) 

GE 


b 

6(6-tl) 

1 

l-\-b — bp 

NB 

rr‘){i7)"{A)‘ 

ba 

b(b + l)a 

(l+6_6p)-“-l 
6a(p—1) 

SF 

l<k<N 

0, k > N 

Tr{l-1) 

Hn 

Lfil) 

see [36] 


TABLE I. The burst-size distribution (BSD), its mean, vari¬ 
ance, and the corresponding f{p) [defined immediately after 
Eq. (4)] for single step influx (SSI), A'-step influx (KSI), ge¬ 
ometric (GE), negative-binomial (NB), and scale-free (SF) 
BSDs. Here, h//'^ is the A'-th harmonic number of order 7 . 


FIG. 2. (Golor online) Panels (a-c): ratio between the MET 
for various BSDs and the MET for the SSI case as a function 
of different characteristic parameters, see Table I. Solid line 
- theoretical result ( 6 ); markers - simulation results. Panel 
(d): the MET for SF distribution as a function of 7 . In all 
panels, we have multiplied our analytical results by constant 
pre-exponential factors obtained by a numerical fit: 21, 14, 10, 
9, 2 for single-step influx, A'-step influx, geometric, negative- 
binomial, and scale-free BSDs, respectively. In all panels ko = 
10 and a = 1.1; panel (a): N = 90; panel (b): N = 105; panel 
(c): A = 105 and a = 2; panel (d): A = 250. 


the reader, we have chosen the BI rate in our model in 
such a way that the deterministic description remains 
unchanged compared to the SSI case. Therefore, this re¬ 
duction in the MET is a net effect of BI and occurs due to 
the right tail of the BSD. The latter encodes occasional 
large bursts of individuals, which allow for a more rapid 
escape from the metastable state. 

Let us now discuss in detail the different BSDs we have 
used. We begin with the simple case of influx occurring 
in bursts of constant size AT > 1, see Table I. While this 
BSD differs from the SSI case only by the mean, the fact 
that influx occurs in bursts of size AT > 1 at a rate ko /K 
suffices to reduce the MET exponentially, see Fig. 2a. As 
expected, at AT = 1 the AT-step case reduces to tssi- 

Next, we consider the case of geometric BSD, which 
occurs, e.g., in cell biology [18-21, 39], and also the more 
general case of negative-binomial BSD, whose mean and 
variance can be independently controlled. In both these 
cases, one observes an exponential reduction in the MET 
compared with the SSI case, see Fig. 2b-l-c. As expected, 
see Fig. 2b-|-c, for 6 —>■ 0 for which the width of the geo¬ 
metric and negative binomial BSDs vanishes, see Table I, 
the MET reduces to the SSI case, since f{p) 1. 

Until now we have considered BSDs with an Oil) vari¬ 
ance. To fully explore the effect of BI we have also stud¬ 
ied the case of bursts sampled from a scale-free distribu¬ 
tion, D{k) ~ k~'^, which plays a key role in many pro¬ 
cesses in nature [40]. Here, for 7 < 3, the BSD variance 
scales with the system size. Thus, we apply a natural 
cutoff for the scale-free BSD, see Table I and [41]. In 
this case we also observe excellent agreement between 
the analytical and numerical results, see Fig. 2d [36, 42]. 
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In the bifurcation limit the results drastically simplify. 
Here one can find an analytical expression for the ratio 
of the accumulated actions, (j), as a function of the first 
two moments of the BSD. Let us define the bifurcation 
parameter e = — AaKo 1 , which guarantees that 

the fixed points are close q 2 — qi ~ e. This limit allows for 
the parabolic approximation of the activation trajectories 
for which p — 1 ~ [43]. Expanding g{p) around p = 

1 , we find f{p) = (l/(fc)) Em=i 

where is the m-th derivative of g{p) evaluated 

at p = 1, which yields combinations of BSD moments. 
Retaining the leading- and subleading-order terms yields 

f{p) ~ l-b 2 B(p-l) ; ^ [var(fc)/(fc) -b (fc) - 1 ], ( 8 ) 

where (fc) and var(fc) are the mean and variance of the 
BSD. Here, the leading order term of /(p) corresponds 
to the SSI result, see Table I, while the subleading order 
term encodes the effect of BI for any arbitrary BSD. 

Defining a rescaled momentum coordinate p = (p — 
l)/e^, /S.S [Eq. ( 6 )] becomes AS” = [< 7 +(p) — 

q- (p)]dp. Using Eq. ( 8 ), in the leading order in e -C 1, the 
integrand becomes q+{p) — q-{p) ~ 4Koe\/l^^2^(l+Tl), 
while Pmax — 1/[2(1 + B)]. Plugging these results into 
the integral for AS' and preforming the integration, we 
find AS in the leading order in the bifurcation limit 
e 1: AS^ ~ AS|gj/(l -b B). Here, using Eq. (7), 
ASgsi — 4ko£:^/ 3 is the result for the SSI case in the 
bifurcation limit. As a result, close to bifurcation, the 
ratio of accumulated actions reads 

/ = ASVAS|s, = (1 + B)-i. (9) 

The validity of the results for AS^ and ASggj can be 
checked a-posteriori. Indeed, the WKB theory formally 
requires that NAS ^ 1 [7]. That is, e has to satisfy 
7 V- 1/3 <g; e 1 for the analytical results to be accurate. 
Also, note that due to the Cauchy-Schwarz inequality, B, 
defined by Eq. ( 8 ) satisfies B > 0. Only in the case of 
SSI, for which (k) = 1 and var(A:) = 0, we have B = 0 
and (ffi = 1. Thus, close to bifurcation, we have analyti¬ 
cally shown that BI of any arbitrary BSD exponentially 
reduces the MET. We have numerically confirmed that 
this result also holds beyond the bifurcation limit. 

Using the numerical simulation results from Eig. 2a-c, 
we compare in Eig. 3a the numerical value of (ji with the 
theoretical bifurcation-limit result, (ffi (9). The numeri¬ 
cal results for three different BSDs exhibit a collapse onto 
the theoretical result, demonstrating a universal scaling 
of the exponential reduction in the MET (compared to 
the SSI case) on B, for any BSD. The fact that in Fig. 3a, 
e = 0.75, indicates that the bifurcation-limit result may 
extend beyond its formal region of validity £ <C 1. Fig¬ 
ure 3b shows [AS” — AS^\/AS - the relative error of the 
accumulated action in the bifurcation limit with respect 
to the general case - as a function of e. As expected, the 



10~^ 10~’ l(f 

e 


FIG. 3. (Color online) Panel (a): The accumulated action 
ratio (j)j obtained from numerical simulations (markers) com¬ 
pared with the theoretical bifurcation result, (/>*' [Eq. (9)] 
(solid line). Parameters and data are the same as in Fig. 2a- 
c, and e ~ 0.75. Panel (b): Relative error of the accumu¬ 
lated action in the bifurcation limit versus the general case, 
[AS — AS'*’|/AS', as a function of e. Here kq = 0.2 for all 
lines; GE: b = 3; NB: fe = 3, a = 5; SF: N = 1000. 


relative error vanishes as £ —>■ 0. That the bifurcation 
limit holds, for various BSDs, as long as £ <C 1, verifies 
our theoretical assumption, and further demonstrates the 
robustness and applicability of the bifurcation result. 

We have calculated the mean escape time (MET) from 
a metastable state in the case where immigration occurs 
in bursts rather than via arrival of single individuals at 
a certain rate. This was done using an arbitrary burst 
size distribution (BSD). We have shown that the effect of 
bursty immigration can exponentially decrease the MET 
from a metastable state. Close to the bifurcation limit, 
when the metastable state and the corresponding unsta¬ 
ble state (which serves as the barrier for escape) are close, 
we have found a simple analytical expression for the de¬ 
crease in the MET as function of the first two moments of 
the BSD. Our theory was demonstrated on several pro¬ 
totypical examples of BSDs including the geometric and 
negative-binomial distributions, appearing in cell biology, 
as well as scale-free distributions, which describe many 
social, biological, and communication networks. 

We anticipate that step-size noise appearing in prob¬ 
lems of similar nature will also strongly influence the 
long-time behavior of stochastic populations. The for¬ 
malism presented here can be extended to two additional 
types of step-size noise: bursty reproduction and death in 
groups. The former appears in the context of evolution¬ 
ary biology [44, 45], whereas the latter has been recently 
observed in nature in populations of Saiga antelopes [46] . 
It would be interesting to study how, e.g., extinction of 
a population is affected by having a fluctuating offspring 
number per birth event, or by having an instantaneous 
removal of an a-priori unknown number of individuals. 
Here, however, the analysis is expected to be more in- 
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volved than the case of bursty influx, due to the strong 
coupling between demographic and step-size noise [47]. 

Finally, the formalism that we have developed here for 
static step-size noise makes a first step in the systematic 
study of step-size noise with arbitrary correlation time. 
Importantly, the latter may confer new insight as to the 
role of other types of non-demographic noise, such as 
extrinsic or environmental noise, see e.g. Refs. [48, 49], 
on stochastic populations dynamics. 

We thank Baruch Meerson for useful discussions. This 
work was supported by Grant No. 300/14 of the Israel 
Science Foundation. 
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